{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://csdms.colorado.edu/wiki/ESPIn2020\"><img style=\"float: center; width: 75%\" src=\"../../media/ESPIn.png\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advection\n",
    "\n",
    "@Auhtor: Benjamin Campforts\n",
    "\n",
    "### Introduction\n",
    "\n",
    "- **What is advection?** Advection is a transport mechanism of a substance or conserved property by a fluid due to the fluid's bulk motion (Wallace and Hobbs, 2006).\n",
    "\n",
    "- **Examples:**  What can be advected: temperature, moisture, aerosols, ...\n",
    "\n",
    "- **Why discretise?** Because not for all solutions the analytical solution is known, in fact, most of the time it is not known. Especially in the field of earth sciences where sets of simultaneous, nonlinear, partial differential equations has to be solved using known analytical methods. To solve such coupled equations, numerical computational methods are required to discretize and solve the PDE's on a lattice. In this notebook, we will only deal with a finite difference scheme, using the simplest form of a truncated Taylor series expansion. Many alternative methods are available to solve PDE's. Nevertheless, many models in geosciences use finite difference methods, because of their conceptual simplicity and ease of coding (Prof. N. Van Lipzig)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The advection equation: some theoretical background###\n",
    "\n",
    "A concepturally very simple PDE is the Advection Equation: \n",
    "\n",
    "$$\\frac{\\partial C}{\\partial t} = -v \\frac{\\partial C}{\\partial x}\n",
    "\\label{eq:1}\\tag{1}$$\n",
    "or \n",
    "$$\\frac{\\partial C}{\\partial t} +v \\frac{\\partial C}{\\partial x}=0\n",
    "\\label{1b}$$\n",
    "\n",
    "where $C$ is the aerosol concentration and $v$ is a constant windspeed at which aerosol concentrations are advected. The equation above is a prototype of an **initial value problem**: The solution is obtained by using the known initial values and marching or advancing in time. The solution of this equation can be obtained directly from the initial conditions:\n",
    "\n",
    "$$ C(x,t) = C(x-vt,0) \\label{eq:2}\\tag{2}$$\n",
    "\n",
    "In the following we will discretize Eq. $\\ref{eq:1}$ using a first order, finite difference upwind scheme. \n",
    "\n",
    "All FDM's can be represented with a numerical stencil:\n",
    "\n",
    "<img src=\"../../media/Stencil1.png\" style=\"width:3in;height:1in\" />\n",
    "\n",
    "For the upwind scheme in particular, we only need three nodes: \n",
    "<img src=\"../../media/Stencil2.png\" style=\"width:3in;height:1.25in\" />\n",
    "\n",
    "**BUT** the nodes you need to consider depend on the sign of the velocity: \n",
    "<img src=\"../../media/UpwindScheme.png\" style=\"width:6.5in;height:3.76in\" />\n",
    "\n",
    "If $v$ is positive, than the direction of the Aerosol flow is from left-to-right along the x domain. To capture such a positive flow, the upwind discretisation scheme uses the info from the neighbouring node in the upwind direction, that is, the node to the left, indicated by j-1 on the sketch. In this case the the Finite Difference Equation (FDE) of Eq. $\\ref{eq:1}$ can be written as:\n",
    "\n",
    "$$\\frac{\\tilde{C}_j^{n+1} - \\tilde{C}_j^n}{\\Delta t} + c\n",
    "\\frac{\\tilde{C}_j^n - \\tilde{C}_{j-1}^n}{\\Delta x} = 0\n",
    "\\label{eq:3}\\tag{3}$$\n",
    "\n",
    "If however, the velocity is negative, the FDE becomes: \n",
    "\n",
    "$$\\frac{\\tilde{C}_j^{n+1} - \\tilde{C}_j^n}{\\Delta t} + c\n",
    "\\frac{\\tilde{C}_{j+1}^n \\tilde{C}_j^n}{\\Delta x} = 0\n",
    "\\label{eq:4}\\tag{4}$$\n",
    "\n",
    "To obtain this equation, we have taken discrete values for $x$ and $t$: $x_j=j\\Delta x$, $t_n = n \\Delta t$. The solution of the finite difference equation is also defined at the discrete points ($j\\Delta x$,$n \\Delta t$): $\\tilde{C}_j^n=\\tilde{C}(j\\Delta x,n \\Delta t)$. Note that we used $C$ to denote the solution of the PDE\n",
    "(continuous) and $\\tilde{C}$ to denote the solution of the finite difference equation (FDE), a discrete solution. In the remainder of the notebook, we will use the symbol $C$ for both the analytical and numerical solution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The advection equation: Hands on ###\n",
    "\n",
    "In the following we will calculate the temperature anomaly of a discrete, sine shaped temperature perturbation and visualize the result.\n",
    "\n",
    "- Carefully scan each line of the script and try to understand the structure\n",
    "- Run the script to generate a wave \n",
    "- This wave can represent temperature over the 1D domain or a concentration of a pollutant\n",
    "- Play with parameters (amplitude/wave length) to understand its meaning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------------------------------------------------------\n",
    "# Code to calculate the temperature anomaly of a discrete, \n",
    "# sine shaped perturbation and visualise the result\n",
    "# --------------------------------------------------------------------\n",
    "                          \n",
    "# --------------------------------------------------------------------\n",
    "# initialisation\n",
    "# --------------------------------------------------------------------\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# initialisation of model constants\n",
    "dx     = 500      # spatial step in m \n",
    "lam    = 5e4      # wavelength\n",
    "Lx     = 25e4     # size of the domain\n",
    "ampl   = 5        # amplitude (K)\n",
    "\n",
    "# initialisation of variables (including scalars, arrays and matrices)\n",
    "x = np.arange(start= 0, stop = Lx , step = dx)\n",
    "# Logical operations\n",
    "b = x <= lam\n",
    "c = x >  lam\n",
    "# print(c)\n",
    "\n",
    "# --------------------------------------------------------------------\n",
    "# manipulations\n",
    "# --------------------------------------------------------------------\n",
    "# get sine wave\n",
    "y    = np.zeros_like(x)\n",
    "y[b] = ampl*np.sin(x[b]*2*np.pi/lam)\n",
    "y[c] = 0\n",
    "# --------------------------------------------------------------------\n",
    "# Visualisation\n",
    "# --------------------------------------------------------------------\n",
    "# plot the graph\n",
    "plt.figure()\n",
    "plt.plot(x/1.e3, y, '--k') \n",
    "plt.ylim((-6, 6))\n",
    "plt.xlim((0, max(x/1e3)))\n",
    "plt.xlabel('x (km)')\n",
    "plt.ylabel('T (K)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets solve the advection equation using the upwind FDE. You can use either for loops or vector operations (numpy arrays) to iterate over the temperature vector in space. Start with the solution for positive velocity (that is, solve Eq. 3) and assume: \n",
    "\n",
    "- Wind speed $v$ = 10 m/s\n",
    "- Total model run time = 3 hours\n",
    "- calculate dt using the CFL criterion: \n",
    "$dt = \\frac{0.9v}{dx}$\n",
    "- dx = 500 m (see previous code block)\n",
    "- Periodic boundary conditions: the solution at the left boundary depends on the solution at the right boundary if velocity is positive and vice versa for negative velocities. \n",
    "- Plot the wave every $10^5$ iterations \n",
    "\n",
    "The structure of the code should look like this: \n",
    "\n",
    "~~~\n",
    "# data, parameter values\n",
    "v    = \n",
    "totT = \n",
    "dt   = \n",
    "print(dt)\n",
    "\n",
    "# Initialise numpy array \n",
    "T_num = \n",
    "T_num[:] = \n",
    "\n",
    "# Numerical solution\n",
    "it = \n",
    "\n",
    "while it*dt<totT:\n",
    "    it +=1\n",
    "    if v>0:\n",
    "        # here comes the solution of the advection equation  \n",
    "    \n",
    "    if it%1e5==0: \n",
    "        plt.plot(x/1.e3, T_num, 'b') \n",
    "        plt.xlim((0, max(x/1e3)))\n",
    "        plt.xlabel('x (km)')\n",
    "        plt.ylabel('T (K)')\n",
    "        plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, set the wind speed $v$ to -10 m/s and calculate the advection equation for negative velocities using Eq. 4. Copy paste your previous solution and add the negative component to it. Do not forget to reset the Temperature array before you start. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you find an analytical solution for this problem? Implement eq. 2. You need just one line of code (HINT use the `np.mod` statement). Plot and compare with the numerical solution. What is going on here. What can we do to resolve this issue? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice your skills: Eyjafjallajokull- Part 2\n",
    "\n",
    "Back to the problem of the Eyjafjallajokull volcano. This time we will solve the advection diffusion equation in 1D: \n",
    "\n",
    "$$\\frac{\\partial C}{\\partial t} = -v \\frac{\\partial C}{\\partial x} + D\\frac{\\delta^2 C}{\\delta x^2} \n",
    "\\label{eq:5}\\tag{5}$$\n",
    "\n",
    "Copy past your script form the notebook where you solved the case for pure diffusion. \n",
    "Now adjust the script by adding advection. You can use the code derived in the previous blocks. \n",
    "\n",
    "You only need one additional model parameter: the advection velocity $v$. \n",
    "- Assume an advection velocity of 10 km/h (that is, the wind speed at which ash aerosols are advected towards the continent)\n",
    "- Calculate the time step by combining the CFL criterium for advection (see above) and diffusion:\n",
    "~~~\n",
    "dt_a = 0.9*v/dx\n",
    "dt_d = dx*dx/D/2.5\n",
    "dt = min(dt_a,dt_d)\n",
    "print('dt is: ' + str(dt) + 'hours')\n",
    "~~~\n",
    "- Same question: after how many hours do we get 5 ppm ash aerosols in Brussels?. \n",
    "\n",
    "Enter the parameter values and initialize model variables: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code up the while loop. \n",
    "The structure should look like this:\n",
    "~~~\n",
    "it =0 \n",
    "while ...\n",
    "    it+=1\n",
    "    \n",
    "    # Advection \n",
    "    if v>=0:\n",
    "        ...       \n",
    "    elif v<0:\n",
    "        ...\n",
    "        \n",
    "    # Boundary conditions\n",
    "    C[0] = ...\n",
    "    C[-1] = ...    \n",
    "        \n",
    "    # Diffusion \n",
    "    ...\n",
    "    \n",
    "    # Source term\n",
    "    C[ind_vol]...\n",
    "    \n",
    "\n",
    "    if it%100==0:        \n",
    "        plt.plot(x,C)\n",
    "        plt.scatter(x[ind_Bru],C[ind_Bru],c='r')\n",
    "        plt.title('Time is: ' + str(it*dt) + ' sec')    \n",
    "        plt.show()\n",
    "            \n",
    "print('Concentration reached after: ' + str(int(it*dt)) + ' hours')\n",
    "print('or : ' + str(int(it*dt/24)) + ' days')\n",
    "~~~\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's more like it! Congratulations, you have reached the end of this notebook!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
